Pricing options with VG model using FFT 



Andrey Itkin 

Rutgers University, New Jersey, 
Department of mathematics 
aitkin@math.rutgers.edu 

Summary 

We discuss various analytic and numerical methods that have been used to get option prices 
within a framework of the VG model. We show that some popular methods, for instance, Carr- 
Madan's FFT method [1] could blow up for certain values of the model parameters even for an 
European vanilla option. Alternative methods - one originally proposed by Lewis, and Black- 
Scholes-wise method are considered that seem to work fine for any value of the VG parameters. 
Test examples are given to demonstrate efficiency of these methods. Convergency of all methods 
is also discussed. 
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1 Introduction 



This paper summarizes some results of work originally initiated by Peter Carr. It supposes to 
investigate various numerical and analytical methods of option pricing using VG model in order 
to find out which algorithm is most efficient. 

Let us first give a brief overview of the VG model. The Variance Gamma (VG for short) 
process was proposed by Madan and Seneta (see [?]) to describe stock price dynamics instead of 
the Brownian motion in the original Black-Scholes model. Two new parameters: 9 skewness and v 
kurtosis are introduced in order to describe asymmetry and fat tails of real life distributions. The 
VG process is defined by evaluating Brownian motion with drift at a random time specified by 
gamma process. In other words, the VG model with parameter vector (a, v, 9) assumes that the 
forward price satisfies the following equation 



]nF t = \nF + X t + ujt, (1) 

where 

X t = 9 lt (l,v) + aW lt{l ^ (2) 

and 7t(l,^) is a Gamma process playing the role of time in this case with unit mean rate and 
density function given by 



7i(l,^) 



(x) 



X" e " 



(3) 



In the Eq. (1) u is chosen to make Ft a martingale. 

The probability density function for the VG process may be written as 



ht(x) 
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or after integration over g 
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7 



where K is the modified Bessel function of the second kind. The characteristic function <f>j t n^(u) 
for the VG process has remarkably simple form 



r°° i 

cj> t {u) = (E™) = / h t (x)e lux dx = j -. (6) 

Jo (l-iOvu+ -a 2 vu 2 )v 

Another derivation of this expression could be obtained when conditioning on time change like 
in Romano- Touzi for stochastic volatility models 
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£ , x(l — iuv) 

_t [°° yi^e'v _t 
(l-iuv) - / — — — dy = (l-iuv) » . (7) 

Jo ^r(£) 



0x t (n) = E[e iltXt ] = E[E[e iuXt | 7t(l,z/)]] = EfEfe^M 1 ' 1 ^^^'')) | 7f (l,z/)]] 
= E[e™ e,7 *( 1 '^ _ i w2 °' 27 ^ 1 '^] = E[e i ( u0+ ^ u2(j2 ) 7 *( 1 ' i ^] 

= (f> Jt (i jU )(u6 + i^u 2 a 2 ) = ^ - iOvu + -a 2 vu 2 ^J . (8) 

Now, to prevent arbitrage, we need Ft be a martingale, and, since Ft is already an independent 
increment process, all we need is 

E[F f ] = F , (9) 

or 

E[F e Xt+ ^} = F 0( j> Xt {-i)e ut = F . (10) 

This tells us that 

hx6 Xt {-i) -£ki(l-0i/ -\o 2 v) 1 / n 1 , \ 

™ V = 1/ V 2 L = _] n l 1 _ du _ _ a 2 u \ U 

t t v \ 2 / v ; 

Note that from the definition of uj above, in order to have a risk neutral measure for VG model, 
its parameters must obey an inequality: 

!>#+£. (12) 

Note that risk neutral parameters 9,v,a do not have to be equal to their statistical counterparts. 
Accordingly, the characteristic function of the xt = log St VG process is 

</>{u) = - °- — . (13) 

1 — iOvu H — c^wu 2 
2 

Statistical parameters of VG distribution may be calculated from the historical data on stock 
prices. In particular we have to find the values of the parameters 9*,v* and a* such that the 
folloiwng expression is maximized: 
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n 

Y[h Tj ( Xj ), (14) 
where h Tj (xj) are given by Eq.5 and xj are observed returns per time Tj, i.e. Zj = log(Sj/Sj-i). 

2 Pricing European option 

The value of European option on a stock when the risk neutral dynamics is given by Eq. (1) is 

/oo 
h T (x- (r - q + u;)T)W(e x )dx, (15) 
-oo 

where T is time until expiration, q is continuous dividend and W(e x ) is payoff function that has 
the following form 

W(e x ) = {S e x - K) + - call, W{e x ) = (K — S e x ) + - put. (16) 
Direct calculation allows us to derive the put-call parity relation identical to Black-Scholes case 

C = S e- qT - Ke~ rT + P. (17) 

There are several methods to price a European option under the VG model. One method uses 
the closed form solution derived in [2] . Although the expression is analytic it requires computation 
of modified Bessel functions, and hence may not be as fast as we would like our pricing model to 
be. Therefore, FFT method has been widely utilized to obtain a more efficient pricer. Few flavors 
of the FFT method has been previously discussed with regard to the VG model. 

First of all the FFT method of Carr and Madan [1], nowadays almost standard in math finance, 
was applied to the VG model to price the European vanilla option since the characteristic function 
of the log-return process has a very simple form given above. Further we intend to show, that 
unfortunately this method blows up at some values of the VG parameters. 

Mike Konikov and Dilip Madan [3] proposed another interesting method based on the definition 
of the VG process as being a time changed Brownian motion, where the time change is assumed 
independent of the Brownian motion. This method was described in detail in [3] while has not 
been implemented yet. 

Also Mike Konikov and I independently implemented a modification of the FFT method - the 
Fractional Fourier Transform, which is described in detail in [4, 5]. This method usually allows 
acceleration of the pricing function by factor 8-10, while for the VG model it still demonstrates 
same problem as the original FFT. 

Below we discuss why the Carr and Madan FFT approach fails for the VG model. We propose 
another method, which originally has been developed in a general form by Lewis [6], that seems 
to be free of such problems. 

3 C arr- Madan 's FFT approach and the VG model 

Let us start with a short description of the Carr-Madan FFT method. It was worked out for 
models where the characteristic function of underlying price process (St) is available. Therefore, 
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the vanilla options can be priced very efficiently using FFT as described in Carr and Madan [1]. 
The characteristic function of the price process is given by 



<l>(u,t) =E(e 



iuXt > 



(18) 



where Xf = log(St). Note that the above representation holds for all models and is not just 
restricted to Levy models where the characteristic functions have a time homogeneity constraint 
that 4>(u,t) = e - * 1 ^"), where ip(u) is the Levy characteristic exponent. 

Once the characteristic function is available, then the vanilla call option can be priced using 
Carr-Madan's FFT formula: 



„— alog(AT) roo 

C(K,T) = / Re 

7T Jo 



-iv log(K) 



uj(v) 



dv, 



(19) 



where 



u(v) 



-tT , 



(a + l)i,T) 



(20) 



a 2 + a — v 2 + i(2a + l)v 

The integral in the first equation can be computed using FFT, and as a result we get call option 
prices for a variety of strikes. For complete details, see Carr &; Madan paper [1]. 
The put option values can just be constructed from Put-Call symmetry. 



), T = 0.02,0 = 0.01 





Figure 1: European option values in VG 
model at T = 0.02yrs, K = 90, a = 0.01 
obtained with FRFT. 



Figure 2: European option values in VG 
model at T = 0.02yrs, K = 90, a = 0.01 
obtained with the adaptive integration. 



Parameter a in Eq. (19) must be positive. Usually a = 3 works well for various models. It is 
important that the denominator in Eq. (20) has only imaginary roots while integration in Eq. (19) 
is provided along real v. Thus, the integrand of Eq. (19) is well-behaved. 

But as it turned out, this is not the case for the VG model. To show this let us consider the 
European call option values obtained by Mike Konikov by computing FFT of the VG characteristic 
function according to Eq. (19). 

In Fig. 1 the results of that test obtained using the FRFT algorithm are given for strike K = 90, 
maturity T = 0.02 yrs and volatility a = 0.01. It is seen that at positive coefficients of skew ~ 2 
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and coefficients of kurtosis v ~ 0.5 the option value has a delta-function-wise pick that doesn't 
seem to be a real option value behavior. In Fig. 2 similar results are obtained using a different 
method of evaluation of the integral in Eq. (19) - an adaptive integration. Eventually, in Fig. 3 
same test was provided using a standard FFT method. The results look quite different that allows 
a guess that something is wrong with FRFT and the adaptive integration. One could also note that 
this test plays with an option with a very short maturity. Therefore, to let us make another test 
with a longer maturity. In Fig. 4-6 the results of the test that uses same integration procedures, 
but for the option with K = 90, T = 1, a = 1, are presented. It is seen that for longer maturities 
FFT also blows up almost at the same region of the model parameters. Moreover, it occurs not 
only at positive value of the skew coefficient but at negative as well. Thus, the problem lies not in 
the numerical method that was used to evaluate the integral in the Eq. (19), but in the integral 
itself. 

Option values lor FFT, K = 90, T = 0.02,0 = 0.01 Option values (or FFT, K = 90, T = 1.00,0= 1.00 



Figure 3: European option values in VG Figure 4: European option values in VG 

model at T = 0.02yrs,K = 90, a = 0.01 model at T = 1.0yrs,K = 90, a = 1.0 

obtained with FFT. obtained with the FFT. 



Now having expression Eq. (6) for the VG characteristic function let us substitute it and Eq. (20) 
into the Eq. (19) that gives 



C(K,T) oc 



-alog(K)— rT poo 



7T 



K < 



-iv log(K) 



[a 2 + a - v 2 + i(2a + l)v] ( 1 



idvu H — a 2 uu 2 
2 



dv, (21) 



where u = v — (a + l)i. At small T close to zero the second term in the denominator of the Eq. (21) 
is close to 1. Therefore at small T the denominator has no real roots. To understand what happens 
at larger maturities, let us put T = 0.8, v = 0.1, a = 3, a = 1 and see how the denominator behaves 
as a function of v and 0. The results of this test obtained with the help of Mathematica package 
are given in Fig 7. 

It is seen that at v = at positive the characteristic function has a singularity. To investigate 
it in more detail, we assume v = and plot the denominator as a function of a and (see Fig. 8). 
As follows from this Figure in the interval < a < 2 there exists a value of that makes 
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Figure 5: European option values in VG 
model at T = 1.0yrs,K = 90, a = 1.0 
obtained with the FRFT. 



Figure 6: European option values in VG 
model at T = 1.0yrs,K = 90, a = 1.0 
obtained with the adaptive integration. 



the integrand in the Eq. (21) singular. This means that singularity of the integrand can not be 
eliminated, and thus the Carr-Madan FFT method can not be used together with the VG model 
for pricing European vanilla options. Using FRFT or adaptive integration that both are slight 
modifications of the FFT, also doesn't help. 

Note that for the VG model the authors of [1] derived condition which keeps the characteristic 
function to be finite, that reads 




100 



Figure 7: Denominator of the Eq. (21) Figure 8: Denominator of the Eq. (21) 

at T = 0.8, v = 0.1, a = 3, a = 1 as a at T = 0.8, v = 0.1, a = 3,v = as a 

function of v and 0. function of a and 0. 

Also as can be seen, for 0, v and a corresponding to the above mentioned tests a becomes 
negative that doesn't allow using this method to price the options in terms of strike. 



Andrey Itkin 



Pricing options with VG model using FFT 



8 



In order to solve these problems one needs to find another way how to regularize the integrand, 
i.e. eliminate doing it in the way as Carr and Madan did it using a regularization factor e~ ak . 



4 Lewis's regularization 

Another approach of how to apply FFT to the pricing of European options was proposed by Alan 
Lewis [6] . Lewis notes that a general integral representation of the European call option value with 
a vanilla payoff is 

/oo 
(e x - K) + q(x,x ,T)dx, (23) 
-oo 

where x = log St is a stock price that under a pricing measure evolves as St = So exp[(r— q)T+Xx, 
r—q is the cost of carry, T is the expiration time for some option, Xt is some Levy process satisfying 
M[exp(iuXT)] = 1, and q is the density of the log-return distribution x. 

The central point of the Lewis's work is to represent the Eq. (23) as a convolution integral and 
then apply a Parseval identity 

f(x)g(x - x)dx = -L [°° e- iuxo f(u)g(u)du, (24) 
, 2ir J_ 00 

where the hat over function denotes its Fourier transform. 

The idea behind this formula is that the Fourier transform of a transition probability density 
for a Levy process to reach Xt = x after the elapse of time t is a well-known characteristic 
function, which plays an important role in mathematical finance. For Levy processes it is <j>t(u) = 
M[exp(iuXt)],u S 9i, and typically has an analytic extension (a generalized Fourier transform) 
u — > z G C, regular in some strip Sx parallel to the real z-axis. 

Now suppose that the generalized Fourier transform of the payoff function w(z) = e tzx (e x — 
K) + dx and the characteristic function <j>t{z) both exist (we will discuss this below). Then from a 
chain of equalities the call option value can be expressed as follows 



C T (x ,K) = e _rT E [(e x — K) + ] = —- — E / e~ izXT w(z)dz 

^ 7r Uifj.— oo 



rT r rifi+oo 

(25) 



e~ rT 

E 

2tt 



ifi+oo 

e - iz[xo+{r - q+u] )T] e -izX T ^ dz 



tfl — OO 



„—rT riii+oo c ~ r T ri/i+oo 



e iz ^ X0+( - r - q+ul ^E[e- lzXT ]w(z)dz = / e~ izY $x T {-z)w{z)dz. 

Here Y = xq + (r — q + u)T, /i = Im z. This is a formal derivation which becomes a valid proof 
if all the integrals in Eq. (25) exist. 

The Fourier transform of the vanilla payoff can be easily found by a direct integration 

/oo jy-iz+l 
e izx, e x _ K ^+ dx = — _ _ Jmz > 1 ^ 
-oo Z %Z 
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Note that if z were real, this regular Fourier transform would not exist. As shown in [7], payoff 
transforms w(z) for typical claims exist and are regular in their own strips S w in the complex 
z-plane, just like characteristic functions. 

Above we denoted the strip where the characteristic function <p(z) is well-behaved as Sx- 
Therefore, <j>{—z) is defined at the conjugate strip S x . Thus, the Eq. (25) is defined at the strip 
Sy = S x P| S w , where it has the form 

Ke~ rT r i ^+ 00 . , dz 
C(S,K,T) = — / e-^cPxA-z)^—, ^eSy, (27) 

^ 7r Jiii— oo Z — IZ 

and k = \og(S/K) + (r - q + u)T. 

The characteristic function of the VG process has been given by the Eq. (13) and is defined in 
the strip (3 — 7 < Im z < f3 + 7, where 



e / 2 e 2 



P = — > 7 = \ — + — + 2(Rez) 2 . (28) 
a V va a 



This condition can be relaxed by assuming in the Eq. (28) Rez = . Accordingly, 4>{—z) is 
defined in the strip 7 — /3 > Im z > —f3 — 7. 
Now let us choose Im z in the form 




jt = kz = < 1 + + t (29) 



Taking into account the Eq. (12) which makes a constrain on the available values of the VG 
parameters, it is easy to see that fi defined in such a way obeys the inequality fx < 7 — /3. On the 
other hand, as also can be easily seen, jjl > 1 at any value of and positive volatilities a, and the 
equality is reached when = 0. It means, that Im z = fx lies in the strip S x as well as in the strip 
S w , i. e. fj, £ Sy- 

Now one more trick with contour integration. The integrand in Eq. (27) is regular throughout 
S x except for simple poles at z = and z = i. The pole at z = has a residue —Ke~ rT i/(2ir), and 
the pole at z = i has a residue Se~ qT i/(2ir) 2 . The analysis of the previous paragraph shows that 
the strip S x is defined by the condition 7 — (3 > Imz > —(3 — 7, where 7 — /3 > 1, and —j3 — 7 < 0. 
Therefore we can move the integration contour to fi\ G (0, 1). Then by the residue theorem, the 
call option value must also equal the integral along Im z = fi\ minus 2iri times the residue at z = i. 
That gives us a first alternative formula 

Kp~ tT f^i+o f ] 7 

C(S,K,T) = Se-^ - I e-^xri-z)^— (30) 



2tt 



l/Ul— OO 



IZ 



For example, with \x\ = 1/2 which is symmetrically located between the two poles, this last 
formula becomes 



C(S, K, T) = Se~ qT - ^\/5Ke- (r+9)T/2 J™ Re e~ iun <5> (-u - £j 



du 



9 1 



(31) 



In other words, if it is valid at Rez = 0, it will be valid for any Rez 
2 This is because <j>r(—i) — e~" T 
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where k = ln(S/K) + (r — q)T, <&(u) = e luulT 4>x T (u) and it is taken into account that the integrand 
is an even function of its real part. The last integral can be rewritten in the form 

Mu)du, Mu) = -^-^ f -« " |) • ( 32 ) 

This can be immediately recognized as a standard inverse Fourier transform, and by derivation 
the integrand is regular everywhere. Indeed, <px T (—u — i/2)\ u= o = (1 — — ^~ t / v \ therefore 
the denominator vanishes if - =6 + ^-. Now using the Eq. (12) one finds that 9 + ^ > 2h + a 2 or 
^j- < — |. Thus, 9 must be negative to turn the denominator to zero. The last equality could be 
also rewritten as 9 + ^_ < Thus, the denominator vanishes if £ < i.e. f must be negative, 
but it is not! Therefore, the characteristic function in Eq. (32) doesn't have singularity at u = 0. 
Thus, a standard FFT or FRFT method can be applied to get the value of the integral. 

In Fig. 9 -10 the results of the European vanilla option pricing with the VG model conducted 
by using this new FFT method are displayed. Two test has been provided with parameters T = 1 
yr, K = 90,a = 0.1 (Fig. 9) and T = 1 yr, K = 90, a = 0.5 (Fig. 10). It is seen that the option 
value surface is regular in both cases. Zero values indicates that region, where the VG constrain 
Eq. (12) is not respected. The higher values of a and O are the lower values of v are required 
to obey this constraint. Therefore, at higher values of v the model is not defined that produces 
irregularity in the graph. This effect is better observable in Fig. 11 that is obtained by rotation of 
the Fig. 10. The above means that the new FFT method can be used with no essential problem. 
A generalization of this method for FRFT is also straightforward. 

In the region of the VG parameters values where an application of the Carr-Madan FFT 
procedure doesn't cause the problem the results of that method are almost identical to what the 
described above method gives. An example of such a comparison is given in Fig. 12 (my NewFFT 
Matlab code vs Mike's FFT code). It is seen that the difference is of the order of 10~ 7 . 



L 



5 Black-Scholes-wise method 

One more method of regularization of the Fourier kernel for the VG model has been proposed by 
Sepp [8] and is also discussed in [9], [10]. The idea is as follows. 

Given characteristic function <px t ( z ) °f the- model M the price of a European option can be 
expressed as 

1 f [™ e - iulnK e iu[lnS+(r - q+ul)T Ux T (u-i) 



nf = - + — / : , < " TV ' du, (33) 

1 2 27r7„ 0O iu<l>x T (-i) V ' 

1 £ r 00 e - iulnK e iu ^ nS+( - r - q+ ^ T Ux T (u) 

2 2 2iry_ 00 iu 

V M = ^ ^T SqU M _ e ~rT KU M ] 

where £ = 1( — 1) for a call(put). Eq. (33) is a generalization of the Black-Scholes option pricing 
formula. Note that <j>x t (0) = 1 by definition, and (j>x±{~ *) i s a function of time to expiry T and 
parameters of the model only. 

Proof: Assume that (j>x(—z) has a strip of regularity < fj, < 1. First we rewrite Eq. (27) as 
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Figure 9: European option values in VG 
model at T = l.Oyr, K = 90, a = 0.1 
obtained with the new FFT method. 



Figure 10: European option values in VG 
model at T = 1.0yrs,K = 90, a = 0.5 
obtained with the new FFT method. 



CdII Op: on Value - VG m 




Figure 11: European option values in VG model at T = l.Oyr, K = 90, a = 0.5 obtained with the new 
FFT method (rotated graph). 



C(S,K,T) 



Ke 



rT fifi+oo 



2ir 



e~ lzk <t> Xl 



dz 



Ifl — OO 



Ke 



-rT 



2vr 



l/l—OO 



*<i>x T {-z) 



z — IZ 
idz 



(34) 



z 



ifi+oo 

l/l—OO 



e-^ k XT (-z) — 



Ke 



-rT 



.(n(h)-n(h)) 



2vr 

In order to evaluate I\ we employ a contour integral over the contour given by 6 parametric 
curves (see Fig. (13): Y\ : z = u,u E (q, R), q, R > 0; : z = R + ib,b £ (0, v); T3 : z = u + iv,u G 
(R, -R); r 4 : z = -R + ib,b £ (v, 0); T 5 : z = u, u £ (-R, -q); F 6 : z = qe ie , 6 £ (vr, 0). As the 
integrand is analytic on this contour we can apply the Cauchy theorem. Also note that the integral 
along curve Tq is a half of the integral along the whole circle around zero which in turn is equal 
to 2iri 2 Res{e~ lzk 4>t{— z) / z) . As the integrals along vertical lines vanish at R — > 00 and at q — > 
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Figure 12: The difference between the European call option values for the VG model obtained with 
Carr-Madan FFT method and the new FFT method. Parameters of the test are: S = 100, T = 
0.5yr, a = 0.2, v = 0.1, = —0.33, r = q = 0. at various strikes). 





Irn z 






V ■ 




V 


-R <f 





q 


R 



Re z 



Figure 13: Integration contour for TZ(Ii) 



the integral along the real axis tends to an integral from — oo to oo, eventually changing variable 
u — > —u we obtain 

K{h) = 7T + I" e -i^K^SHr-^)T]^M du _ (35) 
J-oo 

To compute the 1Z(h) we use a similar contour build around the point z = i, i.e. F± : z = 
u + i, u 6 (q, R), q, R > 0; T 2 : z = R + ib, b G (1, 1 + v); T 3 : z = u + i(l + v), u G (R, -R); T 4 : z = 
— R + ib,b G (v , 1); T$ : z = u + i,u G (-R, — q); : z = i + qe %e , G (0, 7r). Again taking limits 
i? — > oo and g — >• 0, changing variable u — > u — i, we obtain 

K(I 2 ) = ^e^ T (iT + I" e -™lnK e iuVnS+(r- q +u)T] <l>X T (u ~ i) \ 
K V J-oo lU(j)x T {-l) ) 

Substituting these integrals into the Eq. (34) we obtain the Eq. (33) ■. 

The difficulty in using FFT to evaluate the Eqs. (33), as noted by Carr and Madan is the 
divergence of the integrands at u = 0. Specifically, let us develop the characteristic function 
<f>x t (z) with z = u + iv as Taylor series in u 
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(j)X t {z) = E[e~ vXt ] + iuE[xe 



-vXt] 



1 



u 2 E\x 2 e- vXt ] + ... 



(37) 



In Eq. (30) we have to chose z = u — % in the first expression, and z = u in the second one. 
As it is easy to check in both cases that the leading term in the expansion under both integrals 
is l/(iu) which is just a source of the divergence. The source of this divergence is a discontinuity 
of the payoff function at K = St- Accordingly the Fourier transform of the payoff function has 
large high-frequency terms. The Carr-Madan solution is in fact to dampen the weight of the 
high frequencies by multiplying the payoff by an exponential decay function. This will lower the 
importance of the singularity, but at the cost of degradation of the solution accuracy. 

As the Eqs. (33) can be used whenever the characteristic function of the given model is known, 
we can apply it to the Black-Scholes model as well that gives us the Black-Scholes option price 
V B which is a well known analytic expression. Now the idea is to rewrite representation of the 
option price in the Eqs. (33) in the form 



yM = yyM _ yBS} + yBS_ 

The term in braces can now be computed with FFT as 



(38) 



n 



M-BS 



yrM-BS 
L1 2 



v M -V 



BS 



2^ 



2vr 



fe^-i)^ -4>Bs{u-i)e 



vu 



■du, 



(39) 



3 e-™ [<t> Xt {uV™ T ~ <Pbs(u)] 



M-BS 



I'll 



e~ rT KU!f- BS 



du, 



where k = ln(K/S) — (r — q)T, 4>bs(u) = exp y— ^-^-u 2 ^ and 4>x T {—i) = e wT . This is possible 
because we have removed the divergence in the integrals. In addition the magnitude of 4>x t ( z ) ~ 
<Pbs( z ) is smaller than that of 4>x t ( z ) that increases accuracy of the solution. 

In more detail, first terms of the expansion of 4>x t (u)e %uulT — (f>Bs( u ) an d 4>x t { u ~ i)e % ^ u ~^ u ' T 

4>bs(u 



i e 



2 T in series at small u are 



D\\ u =o 

D2\u=0 



(f>x t (u)e 



0x t (u 



- 4> BS {u) = T(9 + u + —)iu + 0(u 2 ) 
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i(u—i)uiT 



4>BS(U 



■\ —S—T 

^)e 2 



a 2 + 



6 + a 2 



-l + u{9 + a 2 /2) 



(40) 

co ) iu + 0{u 2 ) 



However, an usage of these expressions in the Eq. (39) together with the FFT method produces 
an error of the order of 0(u). That is why it is better to choose a small u = e, for instance e = 10~ 6 , 
then computing integrands in the Eq. (39) exactly and substituting Di i 2|«=o ~ D\p\u=e- 

Fig. 14, 15 show the results of our computation of the European option values under the 
VG model. Difference between the Carr-Madan solution and Black-Scholes-wise solution with 
D 1>2 (u = e) and D 1>2 (u = 0) at T = 1.0yr,a = 0.1,6* = 0.1, v = 0.1 are plotted for 200 strikes. It 
is seen that for the first method the difference is of the order of 0.5%. 
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Difference between VG prices. "Carr-Madan" - "BS-wise" 
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Figure 14: European option values in 
VG model. Difference between the Carr- 
Madan solution and Black-Scholes-wise 
solution with D\^,{ u — 0) at T — 
1.0yr,a = 0.1, 6 = 0.1, v = 0.1 



Figure 15: European option values in 
VG model. Difference between the Carr- 
Madan solution and Black-Scholes-wise 
solution with Di^{u = e) at T = 
l.0yr,a = 0.1,0 = 0.1, v = 0.1 



6 Convergency and performance 

Artur Sepp reported in [8] that the convergency of the Black-Scholes-wise method is approximately 
3 times faster than that of the Lewis method. It could be understood because as we mentioned 
above in the limit of small u the difference between the VG solution and the Black-Scholes formula 
which is under the Fourier integral is of the second order in u while in the Lewis method it is of the 
zero order. In other words using the Black-Scholes-wise formula allows us to remove a part of the 
FFT error instead substituting it with the exact analytical solution of the Black-Scholes problem. 

We also fulfilled investigation of how all three methods converge for the VG model. The 
results are given in Fig. 16,17,18. We display log 10 difference between the option price obtained 
with N = 8192, and that with N = 4096, 1024, 512, 256. We don't see much difference in the 
convergency of the Lewis and Black-Scholes-wise method while the Carr-Madan methods behaves 
better at low N. In Fig. 19 we also present the ratio (Cjv=8192 — Cjv=4096)/CiV=8i92 for all three 
methods. The Carr-Madan still converges better for out of the money spot prices while convergency 
of two other methods is similar. 

Cont and Tankov also analyze the Lewis method. They emphasize the fact that the integral in 
the Eq. (30) is much easier to approximate at infinity than that in the Carr-Madan method, because 
the integrand decays exponentially (due to the presence of characteristic function). Bowever, the 
price to pay for this is having to choose m. This choice is a delicate issue because choosing big [i\ 
leads to slower decay rates at infinity and bigger truncation errors and when /ii is close to one, the 
denominator diverges and the discretization error becomes large. For models with exponentially 
decaying tails of Levy measure, [i\ cannot be chosen a priori and must be adjusted depending on 
the model parameters. 

Carr and Madan in [1] compare performance of 3 methods for computing VG prices: VGP 
which is the analytic formula in Madan, Carr, and Chang; VGPS which computes delta and the 
risk-neutral probability of finishing in-the-money by Fourier inversion of the distribution function, 
i.e. according to the Eq. (33); VGFFTC which is a Carr-Madan method using FFT to invert the 
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Figure 16: Convergency of the Black- 
Scholes-wise method. Difference between 
the option price obtained with N = 8192, 
and that with N = 4096, 1024, 512, 256 



Figure 17: Convergency of the Lewis 
method. Difference between the option 
price obtained with N = 8192, and that 
with N = 4096, 1024, 512, 256 



dampened call price; VGFFTTV which uses FFT to invert the modified time value. The results 
are given in Tab. (1). The computation times for the first two methods involve 160 strike levels. 
The first 4 rows of Tab. (1) display 4 combinations of parameter settings, while the last 4 rows 
show computation times in seconds. 





case 1 


case 2 


case 3 


case 4 


a 


.12 


.25 


.12 


.25 


V 


.16 


2.0 


.16 


2.0 


e 


-.33 


-.10 


-.33 


-.10 


T 


1 


1 


.25 


.25 


VGP 


22.41 


24.81 


23.82 


24.74 


VGPS 


288.50 


191.06 


181.62 


197.97 


VGFFTC 


6.09 


6.48 


6.72 


6.52 


VGFFTTV 


11.53 


11.48 


11.57 


11.56 



Table 1: CPU times for VG pricing. Rep- 
resented from [1]. 





case 1 


case 2 


case 3 


case 4 


a 


.12 


.25 


.12 


.25 


V 


.16 


2.0 


.16 


2.0 


e 


-.33 


-.10 


-.33 


-.10 


T 


1 


1 


.25 


.25 


Lewis 


0.031 


0.031 


0.031 


0.031 


Carr-Madan 


0.047 


0.047 


0.032 


0.032 


BS-wise 


0.078 


0.078 


0.062 


0.062 



Table 2: CPU times for VG pricing. Our 
calculations. 



It is seen that the analytic formula is slow while the slowest (and least accurate in case 4) 
method inverts for the delta and for the probability of paying off. 

However, this is not true if one uses a modified method given in the Eq. (39). Our calculations 
show that the performance of the Lewis method is same as the Carr-Madan method, and the 
performance of the Black-Scholes-wise method is only twice worse (because we need 2 FFT to 
compute 2 integrals) (see Tab. 2). 
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Convergency of the Carr-Madan method for VG model 




Figure 18: Convergency of the Carr- 
Madan method. Difference between the 
option price obtained with N = 8192, and 
that with N = 4096, 1024, 512, 256 



Convergency of three methods for VG model 




Figure 19: Convergency of all three meth- 
ods 



7 Conclusion 

We discussed various analytic and numerical methods that have been used to get option prices 
within a framework of VG model. We showed that a popular Carr-Madan's FFT method [1] 
blows up for certain values of the model parameters even for European vanilla option. Alternative 
methods - one originally proposed by Lewis, and Black-Scholes-wise method were considered that 
seem to work fine for any value of the VG parameters. Convergency and accuracy of these methods 
is comparable with that of the Carr-Madan method, thus making them suitable for being used to 
price options with the VG model. 
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